Quien Está familarizado con las posibilidades y funcionalidades del paquete ggplot2, le será fácil usar geometría espacial. Aquí ampliaremos los aspectos relacionados con la representación cartográfica utilizando algunas de esas funciones y otras nuevas. La combinación de la gramática de ggplot2 con la gestión de la información espacial que hace sf, produce un lenguaje de visualización muy intuitivo y escalable. Esta interacción se lleva a cabo principalmente con la función geom_sf(), que es aplicable a todos los tipos de geometrías en el paquete de ggplot2.
sf es compatible con funciones muy conocidas del paquete dplyr (tidyverse).
Para ilustrarlo con un ejemplo sencillo, cargaremos los paquetes tidyverse y sf y leemos un dataset que contiene información sobre los husos horarios. Además, cargamos otros paquetes que usaremos en esta sección.
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1.9000 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
Vistas las ventajas de geom_sf(), usaremos esta función en los siguientes ejemplos. Tenemos a nuestra disposicón el paquete osmdata que nos permite crear y hacer las consultas directamente desde el entorno de R a Open Street Maps. Aún así, el uso de la overpass-turbo.eu puede ser útil cuando no estamos seguros de lo que buscamos o tenemos alguna dificultad en construir la consulta.
Antes de crear una consulta, debemos conocer qué podemos filtrar. La función available_features( ) nos devuelve un listado amplio de las características disponibles en OSM que a su vez tienen diferentes categorías (tags). Están disponibles más detalles en la wiki de OSM aquí. Por ejemplo, la característica shop contiene como categoría entre otros supermarket, fishing, books, etc.
# A tibble: 6 × 2
Key Value
<chr> <chr>
1 shop [[ Too many Data Items entities accessed. | hunting ]]
2 shop agrarian
3 shop alcohol
4 shop anime
5 shop antiques
6 shop appliance
En la primera parte de la consulta debemos indicar el lugar donde queremos extraer la información. La función getbb( ) crea un rectángulo de selección para un lugar dado, buscando el nombre. La función principal es opq( ) que construye la consulta final. Añadimos con la función add_osm_feature( ) nuestros criterios de filtro. Existen varios formatos para obtener el resultado de la consulta. La función osmdata_*( ) envía la consulta al servidor y en función del sufijo * sf/sp/xml nos devuelve el formato simple feature, spatial o XML.
Como ejemplo de puntos de interés (PDIs) de la base de datos de Open Street Maps (OSM), extraeremos las gasolineras.
#rectángulo de selección para la Península Ibéricam <-c(-10, 30, 5, 46)#construcción de la consultaq <- m |>opq(timeout =60*100) |>add_osm_feature("amenity", "fuel")str(q)
Object of class 'osmdata' with:
$bbox : 30,-10,46,5
$overpass_call : The call submitted to the overpass API
$meta : metadata including timestamp and version numbers
$osm_points : 'sf' Simple Features Collection with 70914 points
$osm_lines : 'sf' Simple Features Collection with 101 linestrings
$osm_polygons : 'sf' Simple Features Collection with 7894 polygons
$osm_multilines : NULL
$osm_multipolygons : 'sf' Simple Features Collection with 58 multipolygons
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
bound <-group_by(pi, CNTR_CODE) |>summarise() |>ms_innerlines() plot(bound) # frontera Portugal, España
gs <-st_filter(gas_station$osm_points, pi) # filtramos las estaciones de la PI#mapa final del resultadoggplot(gs)+geom_sf(data = pi, fill ="grey90", color ="white") +geom_sf(colour ="#08519c",alpha = .1,size = .5) +geom_sf(data = bound, linewidth = .8, color ="white") +theme_void()
Ejemplo: Red hidrográfica
Para un ejemplo de un mapa de líneas, representaremos las redes de ríos en la cuenca del Ebro. La siguiente parte es opcional dado que hemos preparado una subset de España a partir de la versión Europea de HydroRivers (V10). Descarga de la base datos original: https://www.hydrosheds.org/products/hydrorivers
# capas disponibles en una gdbst_layers("./data/HydroRIVERS_v10_eu.gdb") # geodatabase # importar la rivers <-st_read("./data/HydroRIVERS_v10_eu.gdb")# límites Españaesp <-esp_get_country(epsg ="4326")# filtramos a Españarivers_esp <-st_crop(rivers, esp)
La la idea es mapear el grosor de línea y la tranparencia según el caudal medio.
# importamos los ríos de Españarivers_esp <-st_read("./data/rivers_espHydrov10.geojson")
Reading layer `rivers_espHydrov10' from data source
`C:\Users\xeo19\Dropbox\RqueR\data\rivers_espHydrov10.geojson'
using driver `GeoJSON'
Simple feature collection with 36807 features and 15 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: -9.472917 ymin: 36.02292 xmax: 4.320511 ymax: 44.08104
Geodetic CRS: WGS 84
# límites de la cuenca de Ebrohyd_ebro <-st_read("./data/cuencas_ebro.geojson")
Reading layer `cuencas_ebro' from data source
`C:\Users\xeo19\Dropbox\RqueR\data\cuencas_ebro.geojson' using driver `GeoJSON'
Simple feature collection with 1 feature and 36 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -4.404174 ymin: 40.33515 xmax: 2.015186 ymax: 43.18007
Geodetic CRS: WGS 84
# filtramos la cuenca rivers_esp <-st_filter(rivers_esp, hyd_ebro) # mapeamos el caudal medio por tamaño de línea y transparenciaggplot(rivers_esp, aes(linewidth = DIS_AV_CMS, # caudal medioalpha = DIS_AV_CMS)) +geom_sf(colour ="#2171b5") +scale_linewidth(range =c(.01, 3)) +# rango de grosorscale_alpha(range =c(.3, 1)) +# rango de transparencialabs(linewidth =NULL, alpha =NULL, title =expression("Caudal medio "(m^2/s))) +theme_void()
Un mapa de símbolos proporcionales usa el tamañó de un símbolo, habitualmente círculos para representar una cantidad. Volvemos a usar el paro registrado por la EPA en el último semestre del año pasado a nivel provincial.
# importar datos de paro EPA INEparo <-read_csv2("data/3996.csv") |>clean_names()
ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
Rows: 5724 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ";"
chr (4): Sexo, Provincias, Tasas, Periodo
dbl (1): Total
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# filtramosparo_prov <-filter(paro, sexo =="Ambos sexos", provincias !="Total Nacional", tasas =="Tasa de paro de la población", periodo =="2023T4") |> dplyr::select(provincias, total) |>mutate(prov_cod =str_extract(provincias, "[0-9]{2}"))# obtenemos las provinciasprov_esp <-esp_get_prov() # unimosparo_prov_sf <-left_join(prov_esp, paro_prov, by =join_by(cpro == prov_cod)) |>st_centroid()
Warning: st_centroid assumes attributes are constant over geometries
ggplot() +geom_sf(data = prov_esp, fill ="grey95", colour ="white") +geom_sf(data = paro_prov_sf, aes(size = total)) +labs(size ="Tasa", title ="Paro en el T4 2023") +theme_void() +theme(legend.position ="top",legend.justification ="left",plot.title =element_text(margin =margin(b =10)))
ggplot() +geom_sf(data = prov_esp, fill ="grey95", colour ="white") +geom_sf(data = paro_prov_sf, aes(size = total)) +scale_size(range =c(.5, 10)) +labs(size ="Tasa", title ="Paro en el T4 2023") +theme_void() +theme(legend.position ="top",legend.justification ="left",plot.title =element_text(margin =margin(b =10)))
ggplot() +geom_sf(data = prov_esp, fill ="grey90", colour ="white") +geom_sf(data = paro_prov_sf, aes(size = total, color = total)) +scale_size(range =c(.5, 10)) +# por áreascale_color_binned_sequential(palette ="SunsetDark") +guides(color=guide_legend()) +labs(size ="Tasa", color ="Tasa", title ="Paro en el T4 2023") +theme_void() +theme(legend.position ="top",legend.justification ="left",plot.title =element_text(margin =margin(b =10)))
Ejemplo: Mapa de accesibilidad media por carretera
En este caso importamos una geometría vectorial con las provincias de España, que contienen datos sobre la accesibilidad a los núcleos de población en 2015. Estos datos representan, en minutos, el promedio y la varianza de los tiempos de recorrido por carretera entre los núcleos de población (> 1500 habitantes por km2) de cada provincia.
prov_acces <-st_read("./data/accesibilidad.shp")
Reading layer `accesibilidad' from data source
`C:\Users\xeo19\Dropbox\RqueR\data\accesibilidad.shp' using driver `ESRI Shapefile'
Simple feature collection with 59 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -18.15839 ymin: 27.63964 xmax: 4.31205 ymax: 43.75056
Geodetic CRS: WGS 84
ggplot(prov_acces) +geom_sf(aes(fill = m))
Para cambiar el color de los límites provinciales a blanco definimos, dentro de la función geom_sf() el parámetro colour = "white".
p <-ggplot(prov_acces) +geom_sf(data = world, fill ="black", colour ="white") +geom_sf(aes(fill = var), colour ="white") +coord_sf(xlim =c(-20, 5), ylim =c(27, 45))p
Aquí hacemos una pequeña paréntesis sobre dibujar límites internos de países. Aplicamos un con conjunto de funciones muy interesante a la ahora de simplificar geometría más.
Y para cambiar la gama de colores de la variable usamos el paquete RColorBrewer, que contiene escalas de color predefinidas ampliamente utilizadas en cartografía (ver colorbrewer2.org). La función brewer.pal permite seleccionar el número de colores y la escala. Ponemos la leyenda en la parte inferior del gráfico con el parámetro legend.position = "bottom", consulta la ayuda de la función theme() y verás que tiene decenas de parámetros para modificar el aspecto del gráfico.
# paleta de coloresp <- p +scale_fill_distiller(palette ="RdYlGn") +labs(fill ="Tiempo (minutos)") +theme_bw() p
En ggplot también existe la posibilidad de modificar la leyenda con la función guides(). En ese caso tendríamos que usar los parametros de aes() (fill, size, colour, shape, etc.) para configurar la leyenda con datos continuos guide_colourbar() o con datos discretos guide_legend(). En la función guide_colourbar() se especifican varios elementos de estilo, desde la posición de las etiquetas hasta la anchura de la barra.
# cambiamos el estilo de la leyendap +guides(fill =guide_colourbar()) +theme(legend.position ="bottom",legend.text.position ="bottom",legend.key.height =unit(.5, "lines"),legend.key.width =unit(3, "lines"),legend.title.position ="top")
Al usar la clase sf no sólo podemos usar la función de geometría geom_sf(), sino también coord_sf(), que nos sirve para definir la proyección o limitar el área de visualización. Es decir, que es posible ajustar la proyección deseada en el proceso de la visualización sin tener que modificar los objetos espaciales originales.
# Cambiamos la proyección a LAEAp +coord_sf(crs =3035, default_crs =4326,xlim =c(-20, 5), ylim =c(27, 45))
Coordinate system already present. Adding new coordinate system, which will
replace the existing one.
Por ejemplo podríamos proyectar, en el proceso de visualización, un mapa global sin la Antarctica a la proyección eurocéntrica LAEA (ETRS89 Lambert Azimuthal Equal-Area).
Como área de interés usamos Suiza en este ejemplo. Con excepción de los límites de lagos, los datos necesarios los obtenemos a través de APIs usando diferentes paquetes. El paquete giscoR permite obtener los límites de países con diferentes resoluciones.
Los límites de los lagos corresponden a una capa de modelos cartográficos digitales (DKM500) que ofrece swisstopo. El objetivo es quedar sólo con los grandes lagos, por tanto excluimos todos aquellos con menos de 50 km2 y también aquellos situados completamente en territorio italiano. Recuerda que con el paquete units podemos indicar unidades y así hacer cálculos.
# importamos los lagossuiz_lakes <-st_read("./data/22_DKM500_GEWAESSER_PLY.shp")
Reading layer `22_DKM500_GEWAESSER_PLY' from data source
`C:\Users\xeo19\Dropbox\RqueR\data\22_DKM500_GEWAESSER_PLY.shp'
using driver `ESRI Shapefile'
Simple feature collection with 596 features and 14 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 2480000 ymin: 1062000 xmax: 2865000 ymax: 1302000
Projected CRS: CH1903+ / LV95
# filtramos a grandes lagossuiz_lakes <-mutate(suiz_lakes, areakm =set_units(SHP_AREA, "m2") |>set_units("km2")) |>filter(areakm >set_units(50, "km2"),!NAMN1 %in%c("Lago di Como / Lario","Lago d'Iseo","Lago di Garda"))plot(suiz_lakes)
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all
Modelo digital de terreno (MDT)
La función get_elev_raster() nos permite descargar un MDT de cualquier región del mundo a través de diferentes proveedores en formato de ráster. Por defecto usa AWS. Un argumento esencial es la resolución que depende de la latitud, la que se puede indicar como nivel de zoom (véase la ayuda). Por ejemplo, aquí usamos nivel 10 que a una latitud de 45º correspondería a aproximadamente 100m.
Después de obtener el MDT de Suiza debemos enmascarar los límites del país. La clase del objeto es RasterLayer del paquete raster, no obstante, el nuevo estándar es terra con la clase SpatRaster. Por eso lo convertimos y después aplicamos la máscara. Finalmente reproyectamos al sistema de coordenadas de Suiza obtenido de los datos vectoriales.
# obtenemos el mdt con mdt <-get_elev_raster(suiz, z =7)
# convertir a terra y enmascarar la zona de interésmdt <-rast(mdt) |>mask(vect(suiz)) # reproyectamos el mdtmdt <-project(mdt, crs(suiz_lakes))# reproyectamos los limite de lagossuiz <-st_transform(suiz, st_crs(suiz_lakes))
Antes de calcular el efecto de sombra, creamos un simple mapa de relieve. En ggplot2 empleamos la geometría geom_raster() indicando la longitud, latitud y la variable para definir el color. Añadimos los límites de los lagos usando geom_sf() dado que se trata de un objeto sf. Aquí sólo indicamos el color de relleno con un azul claro. Con ayuda de scale_fill_hypso_tint_c() aplicamos una gama de colores correspondientes a un relieve, también llamado tintas hipsométricas, y definimos los cortes en la leyenda. En el resto de funciones hacemos ajustes de aspecto en la leyenda y en el estilo del gráfico.
Internamente lo que hace geom_spatraster() es pasar el raster a un data.frame. geom_raster() sirve para una rejilla regular y geom_tile() si fuese irregular.
df <-as.data.frame(mdt, xy = T)ggplot() +geom_raster(data = df,aes(x, y, fill = alt)) +theme_void()
Calcular el hillshade
Recordemos que el efecto hillshade es nada más que añadir una iluminación hipotética con respecto a una posición de una fuente de luz para así ganar profundidad. Las sombras dependen de dos variables, el acimut, el ángulo de la orientación sobra la superficie de una esfera, y la elevación, el ángulo de la altura de la fuente.
La información requerida para simular la iluminación es el modelo digital de terreno. La pendiente y la orientación podemos derivar del MDT usando la función terrain() del paquete terra. La unidad debe ser radianes. Una vez que tenemos todos los datos podemos hacer uso de la función shade() indicando el ángulo (elevación) y la dirección (acimut). El resultado es un ráster con valores entre 0 y 255, lo que nos indica sombras con bajos valores, siendo 0 negro y 255 blanco.
# estimamos la pendientesl <-terrain(mdt, "slope", unit ="radians")plot(sl)
# estimamos la orientaciónasp <-terrain(mdt, "aspect", unit ="radians")plot(asp)
# calculamos el efecto hillshade con 45º de elevación hill_single <-shade(sl, asp, angle =45, direction =300,normalize=TRUE)# resultado final hillshade plot(hill_single, col =grey(1:100/100))
Combinar el relieve y el efecto de sombra
El problema para añadir al mismo tiempo el relieve con su tinta hipsométrica y el efecto hillshade dentro de ggplot2 es que tenemos dos diferentes rellenos para cada capa. La solución consiste en usar la extensión ggnewscale que permite añadir múltiples scales. Primero añadimos con geom_raster() el hillshade, después definimos los tonos grises y antes de añadir la altitud incluimos la función new_scale_fill() para marcar otro relleno diferente. Para lograr el efecto es necesario dar un grado de transparencia a la capa del relieve, en este caso es del 70%. La elección de la dirección es importante, de ahí que debemos tener en cuenta siempre el lugar y el recorrido aparente del sol. (sunearthtools).
<SpatRaster> resampled to 501200 cells.
<SpatRaster> resampled to 501200 cells.
Sombras multidireccionales
Lo que hemos visto es un efecto unidireccional, aunque es lo más habitual, podemos crear un efecto más suave e incluso más realista combinando varias direcciones.
Simplemente mapeamos sobre un vector de varias direcciones al que se aplica la función shade() con una elevación fija. Después convertimos la lista de ráster en un objeto multidimensional de varias capas para reducirlas sumando todas las capas.
# pasamos varias direcciones a shade()hillmulti <-map(c(270, 15, 60, 330), function(dir){ shade(sl, asp, angle =45, direction = dir,normalize=TRUE)} )# creamos un raster multidimensional y lo reducimos sumandohillmulti <-rast(hillmulti) |>sum()# multidireccionalplot(hillmulti, col =grey(1:100/100))
# unidireccionalplot(hill_single, col =grey(1:100/100))
Hacemos lo mismo como antes para visualizar el relieve con sombras multidireccionales.
<SpatRaster> resampled to 501200 cells.
<SpatRaster> resampled to 501200 cells.
La técnica de mezcla de colores es muy útil para obtener resultados notables en el efecto de sombreado. Desde hace poco el paquete ggblend ofrece esta posibilidad. Con el objetivo de combinar varias capas, es necesario insertar los objetos geom_raster() y los scale_fill_*() en una lista seperados por coma. Después le sigue el pipe con la función blend("tipo_de_mezcal") al que le sumamos los otros objetos de ggplot2. En este caso aplicamos la multiplicación como forma de mezcla.
